Exclusion processes: short range correlations induced by adhesion and contact 

interactions 
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We analyze the out-of-equilibrium behavior of exclusion processes where agents interact with 
their nearest neighbors, and we study the short-range correlations which develop because of the 
exclusion and other contact interactions. The form of interactions we focus on, including adhesion 
and contact-preserving interactions, is especially relevant for migration processes of living cells. We 
show the local agent density and nearest-neighbor two-point correlations resulting from simulations 
on two dimensional lattices in the transient regime where agents invade an initially empty space 
from a source and in the stationary regime between a source and a sink. We compare the results 
of simulations with the corresponding quantities derived from the master equation of the exclusion 
processes, and in both cases, we show that, during the invasion of space by agents, a wave of 
correlations travels with velocity v{t) ~ t~^'''^. The relative placement of this wave to the agent 
density front and the time dependence of its height may be used to discriminate between different 
forms of contact interactions or to quantitatively estimate the intensity of interactions. We discuss, 
in the stationary density profile between a full and an empty reservoir of agents, the presence of 
a discontinuity close to the empty reservoir. Then, we develop a method for deriving approximate 
hydrodynamic limits of the processes. From the resulting systems of partial differential equations, 
we recover the self-similar behavior of the agent density and correlations during space invasion. 

PACS numbers: 87.18.Gh, 87.10.Ed, 05.10.-a, 87.10.Hk, 87.17.Aa, 87.18. Hf, 89. 75. Da, 89.75.-k 
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I. INTRODUCTION 

Since their introduction in the 1940s [l|, cellular au- 
tomata have become an essential tool to study collective 
behavior in complex systems starting from the individ- 
ual level in many areas of science: fluid dynamics 
reaction-diffusion problems M_, dynamics of glasses a] , 
epidemiology 0, traffic flow |8l-[l0j . In the recent years, 
there have been many applications in biology, both at the 
intracellular and tissue levels [TTI - [l3j . 

Special cases of cellular automata are the exclusion 
processes vifhich have been successfully applied to study 
motility problems when the concentration of agents is 
such that the geometric hindrance which they impose on 
each other cannot be neglected [13, [iBl • In the exclusion 
processes, each lattice site is occupied by at most one 
agent, so that steric effects (hard-core repulsion) are in- 
corporated from the very beginning and lead to nontrivial 
effects even in absence of other interactions [il[i3. ui- 
terior interactions between agents may be added to study 
intermolecular [18i] . intercellular or interindividual rela- 
tionship 

In practice, these interactions are specified using rules, 
that is, expressions for the time rate of jump of each 
agent to another (empty) lattice site as a function of the 



* lascolani@imnc.in2p3.frj 
t badoual@imnc.in2p3.fr 
t |deroulers@imnc.in2p3.fr| 



present content of each lattice site. All what is needed to 
define the model is the lattice geometry, not necessarily 
regular, and the list of rules. 

Although cellular automata were designed to be effi- 
ciently simulated on a computer, it is helpful to supple- 
ment them with a macroscopic description of the col- 
lective behavior of agents, which often takes the form 
of a partial differential equation (PDE). Even in cases 
where simulations are available (usually repeated a large 
number of times to take into account stochastic noise 
(isl [20j). the PDE may be a compact way to give an 
overview of the large scale behavior of the system, to dis- 
tinguish the universal features of a family of related ex- 
clusion processes in the spirit of a RG-like approach |2]| . 
or to classify them [l^ [22] ■ 

If one wants to study the effects of variations of the 
model parameters, solving PDEs is faster than perform- 
ing stochastic simulations, and PDEs are also useful to 
retrieve analytic results. In cases where simulations of 
the exclusion process is intractable because the number 
of agents in a realistic system is too large (the human 
body contains ~ 10^^ cells, the human brain ~ 3.10^^), 
or because the exclusion process has to be embedded in a 
larger system in a multiscale approach, this macroscopic 
description is essential. 

The proof of existence of a PDE describing an exclu- 
sion process has been the s ubj ect of quite involved mathe- 
matical developments |23H29| . Usually, to get an explicit 
expression for the PDE, one uses a simple approximate 
technique based on the Chapman-Enskog expansion [s^l ■ 
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First, a system of coupled ordinary differential equations 
(ODEs) for the average number of agents in each site of 
the lattice, or equivalently the occupation probability or 
density of each site of the lattice, is derived from the rules 
which define the exclusion process (possibly through the 
use of the so-called master equation [3l|). Because of 
agent interactions and exclusion, each equation gener- 
ally involves joint probabilities or correlation functions, 
like the probability that two nearest-neighbor sites are 
both occupied at the same time. Then, a mean-field-like 
approximation is made to express these correlation func- 
tions as product of site occupation probabilities, as if the 
occupations of two sites were statistically independent. 

Finally, assuming that the occupation probability of 
each site is a regular function of the position when the 
lattice step tends to zero, which amounts to say that the 
typical length scale over which this probability varies is 
much longer than one lattice step, a Taylor expansion 
of this function is substituted into the system of ODEs. 
Truncating the result to lowest non- vanishing order, one 
is left with a PDF for the density of agents as a function 
of space and time. For completeness, let us mention that 
such a derivation of a macroscopic model from a discrete, 
microscopic model can be done in many other settings, 
e.g., to quote a few that were used in biological modeling, 
the cellular Potts model (32l - [33 |. cellular automata on a 
disordered lattice [111 , lattice-gas cellular automata [11] , 
and discrete models with forces (37| . 

In the case of exclusion processes where agents can only 
jump to nearest neighboring sites and where the rules in- 
volve only a short-range interaction between agents [11] , 
the PDE often takes the form of a nonlinear diffusion 
equatiori, and the diffusivity depends on the local den- 
sity [2l-[2l Isi-liH. The interested reader will find ex- 
plicit expressions of the nonlinear diffusivity for a large 
number of such exclusion processes in two recent works 
by Fernando et al. [19] and by Penington et al. [42j . 

This simple, mean-field-like approximation works re- 
markably well when the large-scale behavior of the sys- 
tem is of diffusive nature, i.e., when the nonlinear diffu- 
sion coefficient is positive for all local densities of agents. 
On the contrary, a negative value of the diffusivity is the 
sign that the microscopic dynamics tend to form aggre- 
gates or is subject to demixion. In that case, the average 
occupation number of sites varies on the length scale of 
one or a few lattice steps, hence the hypothesis of regu- 
larity of the density as a function of position when the 
lattice step vanishes is inconsistent, this function cannot 
satisfy a PDE, and the approach breaks down. 

However, even when the diffusivity is always positive, 
the agreement between the density of agents predicted 
from the PDE and the density obtained through an av- 
erage over man y si mulations of the exclusion process is 
not perfect 0,|4l], and this has been proven to be due 
to correlations between the occupations of neighboring 
sites [il]. Usually, the next step beyond the mean- field 
approximation is the so-called pair approximation, when 
one keeps track, at the macroscopic level, of both the 



probability of occupation of each site and of joint proba- 
bilities that two sites are occupied at the same time (here- 
after called two-point correlation functions). This has 
been used, in the case of macroscopically spatially uni- 
form systems, in condensed matter physics [45i] . to study 
random walks El], reaction-diffusion problems [13], epi- 
demic models [48j , ecology [H, HI] , and multicellular sys- 
tems [5l|. Recently, Simpson and Baker extended that 
approach to one-dimensional systems that are macro- 
scopically non uniform, as during an invasion process. 
They found an excellent agreement of the macroscopic 
model with stochastic simulations using two-point cor- 
relation functions of site distances by up to two lattice 
steps [11] or ten lattice steps [1^ . 

In this work, we use a similar approach to study spa- 
tially non-uniform systems in two dimensions, keeping 
the macroscopic model simple by using only nearest- 
neighbor two-point correlation functions, and casting it 
in the form of coupled PDE. While the model is simple, 
its results agree much better with stochastic simulations 
than a PDE for the density alone. Its mathematical ex- 
pression as PDE allows us to analyze the self-similarity 
(or scaling) properties of its solutions in the context of 
invasion of space from a source of cells, as in wound- 
healing-like (5^ or migration assay [55| experiments. We 
show that the scaling properties of correlations may help 
to distinguish between several microscopic mechanisms 
not only in theory, but also in experiments. 

We have in mind applications to living cell migration 
processes, which are essential in a number of biologi- 
cal contexts like development, repair, tumor and can- 
cer progression; therefore, we restrict ourselves to the 
family of exclusion processes where the rate of move- 
ment of one cell depends only on the present contacts 
before moving (which may be preserved or lost) and not 
on the future contacts (contacts with cells which will 
be nearest-neighbor only after the move) — "direction 
then interactions" in the terminology of jToJ. This is 
a realistic setting to study contact interactions (adhe- 
sion (2Q|. [53. f56i. l57i and cell-cell communication phe- 
nomena [43, [5^, [5^, [53), disregarding, e.g., chemotaxis 
or quorum sensing. But the method can be extended to 
more general exclusion processes. A potential applica- 
tion is personalized treatments of invasive tumors such 
as glioma [ll], where computer simulations of a mathe- 
matical model fed by patient-specific parameters will help 
providing the best therapeutic strategy, guide surgical re- 
section, radiotherapy or chemoterapy, and so on. There, 
giving an accurate prediction of the amount of infiltrated 
cells in each part of the brain will be essential. 

This family of exclusion processes is introduced in 
Section [Til In Section [ml we derive the usual mean- 
field macroscopic approximation of them, as well as our 
improved macroscopic models. We compare them to 
stochastic simulations of the exclusion processes. In Sec- 
tion IIVI we go to the continuous space limit and study 
the self-similarity behavior of the solutions. Finally, we 
give a discussion and conclusions in Section |Vl 
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II. THE MODEL 

In our exclusion processes, agents move on a fixed d- 
dimensional lattice, each site of which may contain or 1 
agent. For simplicity, we will apply our results only to the 
bidimensional hexagonal tiling (triangular lattice), but 
they can be extended straightforwardly. In a move, an 
agent can only jump to an empty nearest neighbor site. 
Of course, the reality of biological movement is much 
more complicated than this (for instance, cells deform, 
make protrusions etc.), but the aim is to gain access to 
the macroscopic collective behavior for which we believe 
that a too detailed description may be irrelevant because 
many microscopic degrees of freedom will be "forgotten" 
at large scales and will make simulations and computa- 
tions very difficult if possible at all. 

A. Jump rates and interactions 

The definition of a process is completed with the spec- 
ification of the rate (probability of occurrence per time 
unit) of each jump of an agent, which is assumed not 
to depend explicitly on time or position, but only of the 
content of the lattice sites. 

Let i and j be two lattice sites. We denote V{i) — for 
vicinity — the set of the nearest-neighbors of i on the lat- 
tice, V{i) the numbers of these sites, and v{i) the number 
of sites among them which are occupied (0 < v{i) < V{i), 
V{i) = 1^ = 6 on the hexagonal tiling; v{i) may vary with 
time, but not V{i)). Likewise, Ai{i,j) is the set of lat- 
tice sites which are nearest neighbors of both i and j (but 
distinct from i and j) — M for maintained contact — , 
M{i,j) their numbers (2 on the hexagonal tiling), m(i, j) 
the number of full sites between them, Af{i,j) is the set 
of nearest-neighboring sites of i which are neither j nor a 
nearest-neighbor of j — Af for not-maintained — , N(i,j) 
their number and n(i,j) the number of full sites between 
them, Fig. [H 

For simplicity, we study only some of the processes 
where the rate of any jump, say from site i to site j, de- 
pends only on V{i), m{i,j), n{i,j), M{i,j) and N{i,j), 
but our approach can be extended. We consider the sim- 
plest situation, where there is no other interaction than 
exclusion (hard-core repulsion) , to be the case where each 
agent has the same probability to attempt a jump in a 
given time interval, irrespective of the actual occupancy 
state of the surrounding sites; as a consequence, the rate 
of jumping of an agent is proportional to the inverse of 
the number of possible jumps it can do (number of empty 
surrounding sites). Moreover, all possible jumps in the 
lattice have the same probability to occur, and at large 
scales, the occupation probability obeys a linear diffu- 
sion equation while the mean quadratic distance of each 
agent to its departure point grows linearly with time, 
as if agents would not interact at all and do a simple 
random walk [3, [l^ . In order to facilitate the compari- 
son with results obtained on different lattices, we denote 




FIG. 1. Analytical computations and numerical simulations 
are done in a hexagonal tiling. The distance between two 
nearest neighbor sites is A. Gray filled sites are occupied by 
cells; empty sites are filled in white. The cell marked with a 
white cross symbol attempts to move in the direction of the 
arrow, if the hatched site is empty. The neighbor sites A/" 
and M of the marked cell have black and light gray border 
respectively. Cells in the sites with a black border form non- 
maintained links, which break during the jump of the marked 
cell. Cells in the common neighbor sites between the marked 
cell and the hatched site form maintained links preserved dur- 
ing the jump of the marked cell. 

the rate of jump from site i to site j as Tij /V{i) with, 
by choice, T^j = 1 when there is no other interaction 
than exclusion. The processes we consider (special forms 
of Tij due to interactions) are listed below, along with 
their biological motivation. 

Adhesion model. To study the infiuence of cell-cell 
adhesion on cellular migration, Khain et al. introduced 
a model [11] where 

= (1 - g)'"(^^J')+"(*'^\ (1) 

q € [0, 1] being a constant parameter to quantify the 
strength of adhesion (from 0, no adhesion, like in stan- 
dard Simple Symmetric Exclusion Processes (SSEP) to 1, 
impossible movement). It is assumed that adhesion is 
instantaneously gained or lost, i.e., that the time scale 
of a possible dynamics of adhesion, like the recruitment 
of proteins to build up or strengthen focal adhesions, is 
much shorter than the time scale of migration. 

Gap junctional model. To explain experimental re- 
sults about the influence of gap junction communications 
between cells on the migration of some tumoral astro- 
cytes [ssl . [ssj , Aubert et al. introduced an exclusion pro- 
cess where 

=l-p+(2p-l)min[m(i,j),l]. (2) 
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The parameter p G [0,1], assumed constant and com- 
mon to all cells, allows to interpolate between SSEP for 
p — 1 /2, maximal effect of gap junctions for p = 1, where 
no cell will move unless it has a neighbor cell and will keep 
contact with it, and p — 0, where, to the contrary, no cell 
will move if it is not able to break all existing contacts 
(this case is probably of little relevance to biology, but 
was studied in detail in the context of the glassy dynam- 
ics [6ll - l64| ). Gap junctions are short channels passing 
through two touching cell membranes, which enable the 
passage of small molecules and ions. They are formed 
by two connected hemi- channels, one on each cell mem- 
brane. Different types of gap junctions arise from various 
genie expressions, and a gap junction can be functional 
or closed. The parameter p should be a way to take into 
account such variability. Here again, it is assumed that 
the time scale of establishing effective gap junctions is 
much shorter than the time scale of migration and, ad- 
ditionally, that gap junctions are formed and open with 
any nearest neighbor cells and that their number has no 
influence (i.e., maintaining communication with one cell 
has the same effect as maintaining it with many). 

Linear model. To gain an overview, we finally intro- 
duce a more general model, where 



(3) 



^ 0, /3, and 7 < being constant real numbers. Fig. [T] 
We choose a linear expression for T^.j to keep compu- 
tations simple: the purpose of this model is not to be 
faithful to experiments, but to be illustrative. Choosing 
7 = and /? > 0, one gets a behavior similar to the gap- 
junctional model (the jump rate being now dependent on 
the number of cells with which contact is maintained). 
Choosing /3 = 7 < 0, one gets a behavior similar to the 
adhesion model (for small g, q = 1 and /3 = 7 = — g, the 
behavior should be quantitatively the same as for the 
adhesion model). 



B. Choice of boundary conditions 

The models previously proposed can be studied in dif- 
ferent geometries, that is different boundary conditions 
and different numbers and disposition of cell sources and 
cell sinks, to address various aspects of the exclusion pro- 
cesses such as: relaxation to the equilibrium, steady state 
analysis, approach to the steady state, or perennial out- 
of-equilibrium conditions. One of the geometrical set- 
ups we will focus on is directly inspired by the set-up of 
the cancer cell migration process experiment discussed in 
psl [55I Issj l . The experiments consist in placing an ag- 
gregate of cancer cells (a so-called spheroid) in a Petri 
dish with an agar substrate containing suitable nutri- 
tional needs for the sustainment of the cells. Initially 
piled in the spheroid, the cells slowly exit it and start 
to migrate in the outside region where they avoid over- 
lapping. In the same way, we study the evolution of the 



system starting from a completely empty initial condi- 
tion except for a small central region where all the cells 
are placed. In the proposed geometries, the spheroid is 
represented as a source of cells which can never empty 
out and where no empty sites are allowed at any time. 
When a cell leaves the source to enter the system, the 
free tile in the source, previously occupied by the cell, is 
immediately filled up with a new cell. 

To avoid dealing with infinite lattices, we decided to 
add a sink region. It is an empty reservoir where cells 
are taken away from the system and act as if they were 
driven by a strong apoptotic signal putting them to death 
with an infinite apoptosis rate. Therefore, any cell enter- 
ing the sink is destroyed. When the sink or the source 
extend over a set of contiguous sites that create a closed 
path, they become borders separating the space into in- 
dependent subregions. In this work, for simplicity, we 
consider geometries where there is only one region of in- 
terest, always enclosed between a sink and a source. 

Although cells interact only with their nearest neigh- 
bors and travel during each elementary jump a length 
A equal to the distance between two neighbor cells, the 
boundaries may have infiuence on a long distance be- 
cause of the exclusion rule. Actually, starting from an 
initially empty lattice, except for the source, there will 
be different time regimes: first a non-stationary period 
of time, when the population of cells invades the free 
space and the sink has not yet been reached, then an 
out-of-equilibrium steady state with a constant current 
of cells from the source to the sink (up to stochastic fluc- 
tuations). The first period of time is relevant for the 
migration experiments on Petri dishes. 

Given the set of rules specified in the previous sec- 
tion, we analyze the dynamic evolution toward the steady 
state in two types of geometrical configurations. One on 
a cylindrical surface which has more theoretical advan- 
tages, and one on a planar bounded surface to mimic 
more closely in-vitro experiments of cells migration in 
culture dishes, see Fig. [2] In the cylindrical configura- 
tion, the geometry consists of a two dimensional rect- 
angular space with regular hexagonal tiling having peri- 
odic border condition along one direction and the other 
two sides in contact wi th tw o reservoirs: a source and 
a sink respectively. Fig. 2(a) This geometry is particu- 



larly interesting because it is invariant under translations 
along the direction with periodic boundary conditions. It 
makes possible to describe the properties of the two di- 
mensional system such as density and correlations just 
in terms of one single spatial variable: for example, the 
distance from one reservoir. In the cylindrical geome- 
try, the shape of the reservoirs are fixed and the degrees 
of freedom are the two sizes of the cylindrical surface: 
the length of the circumference and the distance between 
the reservoirs. From numerical simulations, we have seen 
that the results are independent of the length of cylin- 
drical circumference, if this is 10 or more lattice tiles A 
(results not shown) . The distance between the reservoirs 
produces some differences in the steady state, and they 
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FIG. 2. Geometrical disposition of the reservoirs. The black 
sites represent the cell source, the hatched sites are the empty 
reservoir, the empty sites are white and the cells are in gray, 
(a) cylindrical geometry of length L with two small arrows 
connecting neighbor sites showing the two main directions 
of the correlation, (b) radial geometry with two concentric 



circumferences showing the borders of the source and the sink 
with radii Tsrc in white and Tsnk in gray respectively. 



will be discussed in the next section. It is important to 
remark that these considerations hold true when the lat- 



tice is oriented as in Fig. 2(a) as well as when it is rotated 
by 90 degrees. The second geometrical configuration is 
a planar regular hexagonal tiled lattice, hereafter called 
radial geometry. In it, we define an origin O represented 
by a generic tile and two radii from O: Vgrc and Vgnk- The 
tiles with distance from O less than or equal to r^rc func- 
tion as source and the tiles with distance from O bigger 



than Vsnk function as sink. Fig. 2(b) The only constraint 
that exists between the two degrees of freedom, Vsrc and 
fsnk, is Tsrc < fsnk- Nevertheless, numerical results prove 



that for radii larger than two lattice steps A, their spe- 
cific values become weakly infiuential on the dynamical 
evolution of the system (results not shown). Also, in this 
geometrical configuration, the steady state can depend 
on the distance between the reservoirs. 



III. EVOLUTION EQUATIONS ON THE 
LATTICE 

In [60| . the authors analyzed similar systems in the 
mean- field approximation, and in [4^, they commented 
about the discrepancy at the steady state between the 
numerical and the analytical results for particular values 
of their interaction parameter close to the sink. On the 
other hand, here, we address the problem beyond the 
mean-field approximation to investigate the behaviors of 
the correlation at short distance among cells defined as: 



(4) 



where i and j are two nearest neighbor sites, and the an- 
gular brackets stand for the average over all possible con- 
figurations. It is important to stress that the previously 
stated rules and the results shown in the next sections 
will hold true not just as a consequence of the particular 
choice of the hexagonal tiling. Indeed, the phenomenon 
of a correlation wave with similar behaviors is a much 
more general result, and it will be present in other kinds 
of lattices, for example: triangular, square, and random 
tilings, under the condition of a more general definition 
of nearest neighbors like, for example, "a tile is a nearest 
neighbor of another if they have a common vertex" . 



A. Master equation and evolution equations for 
the correlation functions — the general case 

Let us define a lattice as a partition of a n-dimensional 
space C M" in z non-overlapping subsets, each of them 
representing a tile of the lattice identified by an index 
i G N. Let iji be the number of cells in lattice site i; 
rji G {0, 1} so that r]f = r/i, and the total number of cells 
in the system at time t is: Z = X^i^i"'?*- The generic 
configuration of the occupancy states is given by vector 
r] = {rii,rj2, ■ ■ ■ ,r]z) of all the numbers rji. We denote 
P(j7, t) the probability that the the process is in the con- 
figuration rj at time t. 

If Wij is the operator that permutes the contents of 
sites i and j: 

Wij : (?7i,...,77j...,?7j,...,?7^) (771, 77^, ry,, ry^), 

we can express the evolution equation for Pit], t) (the 
master equation) as 



1 



P {W,,,rj,t^ - r;,)T,,,(r,)F(r,,t)] (5) 
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where the first term inside the brackets represents the 
flux of probability from the configurations which con- 
tribute to fill the site i, so that 77^(1 — T]i)Tj_i (Wi_jr]j dt 

is the chance of a cell to move from j to i in the in- 
finitesimal time step dt provided that the system is in 
the state Wij-q. Similarly, the second term represents 
the flux of probability to the configurations where the 
site i is empty. This expression is valid for any exclusion 
process; we assume from now on that Ti_j{ri) vanishes if 
i and j are not nearest neighbors on the lattice (i.e., if 

We define the average of generic quantity A{'q) at time 
i as a sum over all possible configurations of the process: 



(A)=^A(,7)i^(,7,i)- 



(6) 



after a jump of a cell have both sites i and j occupied. 
The equations for the two point connected correlation 
function immediately follow from Eqs. (151 1^ [TT|): 

dtC2{i,J,t) = - J- / ^ [Tfc,.(r,)77fe(l - 77,)+ 

\k£V(i) 

- T,,k{vH{l - Vk)]{V3 - P{j,t))) + I ^ 3- (12) 



In the same way, we can express the evolution equation 
of connected correlation functions for any n. 



B. Expressions for our models 



In particular, we are interested in the average multi-point 
density, or correlation function, on n distinct lattice sites 
Zi, /2, . . . , 



Pn{h 



,ln,t) = iVhVh 



(7) 



— the density (ry^) of cells on a single lattice site i will 
simply be written p{i,t) — and in the connected corre- 
lation function: 

Cn{ll,.-.,ln,t) = {[r]ii - p{li,t)][T]i2 - p{h,t)]... 

[vu-p{um (8) 

which vanishes if the occupation numbers on the sites are 
statistically independent. 

Inserting the master equation Eq. ([5]) into the expres- 
sion of the time derivative of p{i,t) yields the evolution 
equation: 



dtp{i,t) = 



1 



V{^) 



(9) 



(10) 



Using the general definition in Eq. ([7]) for n = 2, the 
equation for the time evolution of the two-point density 
function at the two generic sites i and j is: 



\fcev(i) 



T^,kiv)V^i^-Vk)]VJ) + (11) 



where i -O- j is equal to the first term on the right-hand 
side of the equation above, with the roles of i and j ex- 
changed. In Eq. (|lip . the constraints k ^ j are added 
to ensure that all configurations included in the counting 



Let us give the expression of Tij{T]) for the models 
we introduced in Sec. (jni), on the hexagonal tiling. To 
explicitly compute the averaged quantities of the previ- 
ous subsection, we need to express m{i,j) and n{i,j) in 
terms of rj: 



m{i,j) 



J2 '^fc' 

k£M{i,j) 



n{i,j) 



H Vk- 

keAf{i,j) 



(13) 



Consequently, the transition rate for the adhesion 
model in Eq. ([1]) becomes: 



nk+ni+Vr+'iis+riu 



where k and I are the common neighbors of i and j, and 
r, s and w are the neighbors of i which are not neighbors 
of the future position j. In the gap junctional model, 
when intercellular communication through gap junctions 
drives the system dynamics, Eq. ([2]) can be rewritten as: 

Tijiv) =p{vk + m - mm) + (i -p)(i - fyfe)(i - m), 

where the sites / and k identify the two common neigh- 
bors of i and j. In this way, T^j- will be equal to p if 
both, or just one of the sites k and I are occupied by 
other cells, and it will be equal to 1 — p, if both k and I 
are empty. When rji = 1, the cell in i will share gap junc- 
tional links with all the nearest neighbor sites occupied 
by other cells, but only the gap junctions in the direction 
of the site k or I will be maintained functional during the 
cell transition to the site j. Finally, in the linear model 
the transition rate Eq. ^ is: 

Tij{ri) =a + I3{rik + r;;) + 7(??r + + Vw), 

where the indices k,l,r, s and w have the same meaning 
as in the adhesion model. 

For the gap junctional model and for the linear model 
with 7 = 0, the transition rate is invariant under the per- 
mutation of indices i and j: Tij{rf) = TjAv)- For the 
adhesion model, it is not invariant because the sites be- 
longing to Af(i,j) in the expression for Ti,j(ri) are differ- 
ent from the sites belonging to Af(j,i) in TjAv)- When 
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Tijirf) is symmetric with respect to the indices i and 
J, the probabihty distribution of the configurations at 
equihbrium is the uniform distribution, so that, at equi- 
Ubrium, the occupation of one site is statisticaUy inde- 
pendent from the others. 

In the presence of sources and sinks, the general evolu- 
tion equations for the densities and the correlation given 
in Eqs. ^ [Til [HI) no longer hold. Indeed, if one or more 
points of the multi(single)-point density functions are at 
one lattice step from, or belong to a reservoir, it is neces- 
sary to take into consideration that some changes in the 
configurations of the positions of cells from, or toward 
the reservoirs are impossible and must be excluded. For 
example, when a multi(single)-point density function has 
a point at distance A from a sink, no fiux of cells coming 
from the sink contributes to the final configuration. On 
the other hand, if the reservoir is a source, no cell can 
enter it and any flux of cells toward the source must be 
zero. 



From Eq. ([8]), it follows that if one of the r]i is con- 
stant, then the correlations including the site i are zero; 
consequently, it is easy to express the border conditions 
in terms of the multi-point correlation functions with one 
or more of its points belonging to a reservoir. It is im- 
portant to stress that, despite the presence of reservoirs 
that invalidates the general forms of Eqs. (jHl [TI] [T^ . 
substituting the values of the border condition into the 
evolution equations of the density and the correlations 
explicitly derived from Eqs. (O [H]) produces the same 
correct results as if they were obtained from the evolu- 
tion equations with the ulterior constraints due to the 
presence of reservoirs. 

Let us consider the case of the linear model and ex- 
plicitly compute the equations for the density and the 
two-point connected function in the hexagonal lattice. 
Using Eqs. (|9l[T2|). and dropping the explicit dependence 
of the densities on time and unnecessary indices, the sys- 
tem of equations, in the region of interest far from the 
reservoirs, is: 



k<£V{i) 



P(J, k, s) - p{j)p{k, s) - p{j, i, s) + p{j)p{i, s) 



-7 E 

seAf{i,k) 



si£M(i,k) 

pU, i, k, s) - p{j)p{i, k, s) - p{j, i, s) + p{j)p{i, s) 



-7 E 

se7V(fe,i) 



pU, k, i, s) - p{j)p{k, i, s) - p(j, k, s) + p{j)p{k, s) 



(14) 



p{i,j,k)- p{i)p{j,k) = C2{i,j)p{k) + C2{i,k)p{j) + C3{t,j,k), (15) 

p{i,j, k, I) - p{i)p{j, k, I) = C2ii, i)p{k)p{l) + C2{i, k)p{])p{l) + C2(i, l)p{j)p{k) + 

+ C^{i, J, k)p{l) + C3(Z, J, l)p{k) + C3(i, k, l)p{j) + C4{i, J, k, I), (16) 



where the Eqs. (jl51ll6|) show how to rewrite the quantity 
p{i,j,k) and p{i,j,k,l) for the generic sites i,j,k and I 
in terms of connected correlation functions. 



time t, are: 

p{li) = 1, for li in the source 

p{h) = 0, for li in the sink (17) 

Cnih, ■ ■ ■ ,ln) = 0, if any li is in a reservoir. 



The same set of border conditions Eq. ((T7)) can be applied 
The border conditions for the system Eq. (jl4l) , at any in both cylindrical and radial geometrical dispositions of 
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source and sink, and in addition, they hold true for aU 
three models proposed. 

For a lattice with z tiles, in the system of Eqs. ([9l 
I12p . there are z differential equations for the density and 
even more for the two-point connected correlation func- 
tion. If we consider only the correlations between near- 
est neighbor sites, the process is described by a total of 
(1 -|- (i)z ordinary differential equations plus the border 
conditions, where d x z equations describe the time evo- 
lution of the two-point connected correlations in the d 
main lattice directions (for the hexagonal tiling d = 3). 
In the cylindrical geometry, due to the invariance under 
rotations and reflections along the axial direction, for a 
site there are only two independent equations for the cor- 
relations between nearest neighbor sites: Cy and C_i, see 
Fig. 2(a) therefore, the number of equations reduces to 
4L, where L is the minimum number of sites one has to 
travel through to go from one end to the other of the 
cylinder. 



C. Closure approximation 

To solve the set of equations Eqs. p3 [T71) for the lin- 
ear model, it is necessary to express C-i(i^ j, k, t) in terms 
of known quantities. In more general cases such as the 
gap junctional and the adhesion models, finding a solu- 
tion of the system of equations Eqs. ^ [T^l [H]) requires 
to already know all the correlations C„(/i, . . . ,ln) with 
n > 2, or to add other equations to the initial system 
which allows one to determine the unknown quantities. 
One possible way is to express the multi-point correla- 
tion functions with the highest number of points with an 
approximate expression involving only multi-point corre- 
lation functions with less points. This approach is called 
closure approximation. In comparison to the approach 
in articles (43l . |60| , where mean- field approximation was 
adopted and all the correlations were completely disre- 
garded, we take into consideration the short range two- 
point connected correlation functions with the purpose to 
obtain more information about the dynamical evolution 
of the system and to improve the agreement of the ana- 
lytical results with the stochastic simulations. Our clo- 
sure approximations is the following. All the connected 
correlations defined on more than two points are set to 
zero: 



/ri>2 



(18) 



where, as above, /i, ^2, ■ • ■ , are distinct lattice sites, and 
all the two-point connected correlation functions between 
cells at a distance bigger than one lattice step A are set 
to zero: 

{{V^ - {V^)) (%■ - iV,))) = 0, if .7 i V(*),j ^ (19) 

Therefore, only the information relative to the nearest 
cell couples remains, as if all the clusters and structures 



with more than two points would appear in a completely 
random way (conditioned to the values of the local den- 
sity of cells and of the two-point nearest neighbor corre- 
lation functions). 

These approximations are suggested by the rules of our 
exclusion processes, in which the movement of one cell 
is directly influenced only by the presence or absence of 
nearest neighbors. As we shall see, stochastic simulations 
show that the approximations on the system Eqs. ((9l [T2)) 
produced by Eq. ^TE\i are reasonably good for the gap 
junctional model, the linear model when 7 is small, and 
the adhesion model when q is small (weak adhesive inter- 
actions). For instance, simulations of the gap junctional 
model show that, except cell couples moving together, 
no particular structure or big cluster appears. Never- 
theless, neglecting correlations between sites at distances 
of two lattice steps and above is not so unquestionable 
because of both "repulsive" interactions (exclusion) and 
"attractive" interactions (adhesion). On one hand, it is 
known that, in exclusion processes with no other inter- 
action than mere exclusion, there are long range correla- 
tions [1^ |6^, . From stochastic simulations it seems 
that correlations between cells at two and three lattice 
steps apart are much smaller, but not completely neg- 
ligible (results not shown). In fact, the time evolution 
of these correlations and their behaviors resemble those 
of the connected correlation functions between nearest 
neighbors. On the other hand, our models with large val- 
ues of the adhesion parameter (q or 7) exhibit large-scale 
structures of cells. For instance, the adhesion model can 
be mapped onto the Ising model [l^l , and it can be shown 
that a spontaneous phase separation with clustering of 
cells happens at larges values of q. In such situations, 
the success of the analytical approach in reproducing ex- 
actly results from the stochastic simulations is lost. But, 
at least when aggregation is absent or weak, it is still 
possible to retrieve, from the analytical results for the 
connected correlation functions, important information 
regarding the system evolution which is characteristic of 
each type of interaction (see Sec. lIIIDj and Sec. pVI) ). 

The closure approximations, Eqs. (|18[ [T^ can be 
shortly expressed together with the set of border con- 
ditions Eq. (fT7| : 



Pik) = 1, 
p{k) = 0, 

Cn(W, • • ■ J ^Ji) = 0, 



for li in the source 
for li in the sink 
V in a reservoir or V n > 2 
V/i in a reservoir or li ^ V{l2)- 

(20) 

Other kinds of closure approximations can be applied 
in place of Eq. (|18p . One can systematically extend this 
approximation by dealing with correlations at larger dis- 
tances or with more than two sites [5ll-l53^. One can 
also choose a different scheme, as the Kirkwood Super- 
position Approximation used in [H, H^. However, this 
results in dealing with unbounded connected correlation 
functions which diverge when all sites are empty or when 
all are full, and do not satisfy the geometrical and initial 
conditions proposed here. 
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D. Stochastic simulations and discrete equations 

The stochastic simulations consist of a cellular automa- 
ton where cells evolve through a number of discrete time 
steps moving on a hexagonal tiling and performing a se- 
ries of interaction depending moves described by the rules 
in Sec. In the framework of cellular automata, the 
rules and the change of the state of the cells are intended 
to be applied in parallel (synchronously) [g^l; neverthe- 
less, in the proposed exclusion processes, the parallel up- 
date scheme entails the problem of two cells jumping at 
the same time on the same empty site. To avoid any 
ambiguities in the stochastic simulations, the positions 
of the cells are asynchronously updated following a ran- 
dom order time scheme at each time step all the 
cells are chosen once in a new random order and up- 
dated. At equilibrium, for the gap junctional model and 
linear model with 7 = 0, the symmetry of Tij = Tj^i re- 
sults in no correlations, see Sec. lIIIBl On the other hand, 
simulations with cyclic update scheme produce spurious 
correlations at eq uilibrium. Therefore, the random or- 
dered scheme [69| and the random independent scheme 
[Z^llzil are more appropriate for the simulations. Also 
the clocked random waiting time scheme [72;] should not 
introduce spurious correlations, if the waiting time dis- 
tribution has the first two moments finite. 

The update of a generic cell in the site i consists in 
choosing the new site j at random with equal probability 
between all the nearest neighbors V{i)^ and then, if the 
site J is empty, moving the cell in the new position with 
a probability Q — Tij{ri)At. To be well defined, the 
probability Q that a cell jumps during the time interval 
At requires to fix the simulation time step consistently. 
Therefore, we choose: 



At = 



max(I^,,,(77')) 



adhesion model 

gap junctional model 

linear model 



(21) 



where, in the linear model on the hexagonal tiling, the 
normalization factor max^/(Tij(j7')) = a + max(2/3,0) 
and the parameters must satisfy the inequality a > 
-37-min(2/3,0). 

To reduce the stochastic noise, we averaged the out- 
comes of a series of independent stochastic simulations 
with the same initial conditions for each model and each 
parameter relative to the interaction and geometry. The 
system of ODEs, Eqs. (0 [HI [201) with the initial condi- 
tions corresponding to the CA simulations: 



p{U,t = 0) = 0, 

C2{ll,l2,t = 0)^Q, 



for li not in the source 



are numerically integrated using a fourth order Runge- 
Kutta method and compared with the results of the CA. 

Let us show and discuss in detail the results for the 
more interesting cases. 



Gap junctional model — density profile. The 

density profile at different times on a cylindrical geome- 
try for the gap junctional model with interaction param- 
eter p — 1 and p — 0.9 are shown in Figs. 3(a)[ 4(a) 
Cells, which are initially all positioned in the source at 
a; = 0, migrate away from the full reservoir; as time ad- 
vances, some of them move into the empty region on the 
right resulting in the advancing of the front of the density 
profile which gets closer to the sink a.t x = L. The effect 
of the sink is perceived only at large times, when cells ar- 
rive at the empty reservoir. Before that time, the slope 
of the density profile is strictly negative. The solution 
of the ODEs for the density profile is in good agreement 
with the simulations; nevertheless, some differences are 
noticeable on the right part of the density profile where 
p{x, t) <^ 1 due to the closure approximations. 

Gap junctional model — discontinuity. At large 
times, in the position at one lattice step away from the 
sink, the density presents a discontinuity when p = 1. 
This is due to couples approaching the empty reservoirs. 
When they are in contact with it, the next favorable jump 
in the direction of the sink annihilates one cell of the 
couple, leaving the other in stall without functional gap 
junctional connections. It results in a small accumulation 
of cells on the sites just before the empty reservoir. This 
phenomenon is quite evident in the stochastic simulations 
where the slope of the density profile between the sites 
L — 2 and L — 1 becomes positive. From the solution for 
the steady state equations of the gap junctional model, 
the discontinuity also appears when the source and the 
sinks are few lattice steps apart, for example L — 36, 
but there is no inversion of the slope of the density pro- 
file. To obtain such results, it is necessary to increase 
the distance between the reservoirs. In fact, the density 
at a specific distance from the sink decreases when the 
cylinder length L increases and the decrease is slower at 
position L — 1 than at other distances; eventually, for L 
large enough, it wiU hold true that p{L — 2) < p{L — 1), 
Fig. \5\ Hence, even though the strong effect of the clo- 
sure approximations requires an increase of the distance 
between source and sink, taking into consideration the 
nearest neighbors two-point connected correlation func- 
tion in the equations for the evolution of the system is 
enough to analytically reproduce the discontinuity when 
p — 1. As p decreases, the discontinuity in the density 
profile at the last point close to the sink becomes quickly 
less evident and at p — 0.99, it disappears. The rea- 
son is that, in these cases, single cells can jump and do 
not accumulate; therefore, the slope of the density profile 
remains negative at all positions for all times. Fig. 4(a) 

Gap junctional model — correlation. In 

Figs. [3(b)J [4(b)] the correlation obtained with the cellu- 
lar automata is very different from the correlation from 
the solution of the ODEs. The effect of the approxima- 
tions makes it impossible to recover the exact correlation 
values of the stochastic simulations, but one can see that 
the correlations from the analytical model and simula- 
tions share the same properties. The particular shape 
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of the connected correlation function at a given time t, 
Figs. [3(b)|[4(b)| divides the fi gure in three regions. In the 
first region that goes from the source of the cells to the 
point where the curve crosses the x axis, the correlation 
has a negative tail, and this is due to the exclusion pro- 
cess which forbids several cells from occupying the same 
site. This introduces in the system a short range repul- 
sion between cells, especially close to the source where 
cells are crowded resulting in a negative value of the cor- 
relation. The middle region goes from where the corre- 
lation becomes strictly positive to the point where the 
stochastically simulated correlation is indistinguishable 
from zero (meaning that there the error bar of the corre- 
lation cross the x axis). It corresponds to the zone where 
the density starts to become low, and cells begin to feel 
the lack of neighbor cells. For any p > 0.5, cells tend 
to maintain gap junctions with the neighborhood during 
their motion due to the binding interaction, but the more 
the density approaches an approximative value 1/3, the 
harder it is to preserve a contact with other cells. As a 
result, cells feel the crowding effect less, and the repulsive 
effect of the exclusion process is surpassed by the bind- 
ing interaction. It is in this less dense region that, at the 
microscopical scale, it is possible to see couples moving 
together with a tendency not to separate until another 
neighbor gets close enough to form a new gap junction. 
At p = 1, the effect becomes very strong, and cells moves 
only if they maintain at least one functional gap junction 
and as soon as the density drops down, a population of 
moving gap junctional communicating couples and single 
stalled cells appears. The third region goes from the end 
of the previous region to the sink; here, the correlation is 
almost zero and clearly, the sink is yet too far from the 
population of cells to be perceived. 

Let us consider the time evolution of the connected 
correlation. At the beginning, the repulsion between the 
cells is stronger, but with the advancing of time, it de- 
crease. The peak of the correlation moves toward the 
sink with a velocity v{t) = t^^^^ and its height slowly 
decreases. Eventually, the correlation front will reach the 
sink and its peak slowly disappears leaving a very long 
and small negative correlation tail. At very large time, 
when the correlation is at the steady state, the connected 
correlation between two nearest neighbor sites parallel to 
the bases of the cylinder shows some abrupt changes in 
sign in the proximity of the empty reservoirs. Indeed, the 
large negative values of Cy (L—l) and C± {L — 2) indicate 
a small number of cell couples in perfect accordance with 
the explanation of accumulation of single cells just before 
the sink. 

Linear modeL The linear model with 7 = does not 
differ qualitatively from the gap junctional model. The 
former has a smaller correlation than the latter, and for 
a = at the steady state, the density presents a similar 
discontinuity near the sink. This model has the best 
agreement between the simulations and ODEs solutions. 

Adhesion modeL For the adhesion model, the great- 
est differences between the values of the correlation ob- 
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FIG. 3. Profile of |(a)| the density and |(b)| the connected cor- 
relation at different times for the gap junctional model on a 
cylindrical geometry with p=l and L=36. The dots repre- 
sent the results of the stochastic simulations, and the error 
bars are smaller than the size of the symbols. The lines refer 
to the results of Eqs. ((9l [121 120)) . The insets show the scal- 
ing behaviors. In (b) the insets on the left-hand side refer 
to the analytical results, and those on the right are the re- 
spective results obtained from cellular automata simulations. 
The scaling behaviors of the peak and on the negative tail are 
shown on the top row and bottom row of insets respectively. 



tained from the CA simulations and from the solutions of 
the ODEs are in the region close to the source (see Fig. [6]), 
and they increase with the increasing of the parameter 
q. The repulsion produced by the exclusion process is 
almost inexistent in the case of strong adhesion between 
cells, and the correlation is peaked at the position where 
the density is almost 1/2 showing that, at the front of 
the density profile, cells are almost free, and they do not 
form any structures, but diffuse away. In contrast, at 
higher density, cells gather together, and because in part 
of the crowding and of the high number of links, they 
get trapped, which produces strong correlations between 
cells. The correlation produced by the adhesion has a 
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FIG. 4. Profile of |(a)| the density and |(b)| tlie connected 
correlation at different times for the gap junctional model on 
a cylindrical geometry with p=0.9 and L=36. The dots rep- 
resent the results of the stochastic simulations, and the error 
bars are smaller than the size of the symbols. The lines refer 
to the results of Eqs. (|9] 1121 [20l) . The insets show the scaling 
behaviors. In (b) the inset on the left-hand side refers to the 
analytical results, and that on the right shows the respective 
results obtained from cellular automata simulations. 



peaked shape which moves in the direction of the sink 
with a velocity v(t) cx t^^^^ . This is similar to what hap- 
pens in the gap junctional model. On the other hand, its 
behavior is very different in the other two models. Very 
remarkable facts are that the height of the peak is con- 
stant in time, that it is much stronger than in the gap 
junctional case, and that, for large values of q, there is no 
negative correlation tail close to the source region from 
the beginning of the simulation on. The same constant 
peak of the connected correlation can be obtained in the 
linear model with 7 < and (3 < 0. 

The main differences between CA simulations and an- 
alytical models are due to the closure approximations 
which neglect the long range correlations. Observing the 
ratio between the values from the ODEs and from the 
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FIG. 5. Gap junctional model on a cylindrical geometry with 
p=l. The straight line shows p{L — 1) and the dashed line 
represents p{L — 2). From a numerical fits p{L — 1) goes 



approximately as 



and p{L — 2) goes approximately as 



1 



stochastic simulations of the two-point connected corre- 
lation function at different times, which is approximately 
proportional to ln(i) and (depending on the contact 
interaction) , we conclude that such a discrepancy is pro- 
duced by the exclusion rule, a common factor between 
different models. 



IV. THE HYDRODYNAMIC LIMIT 

In this section, we show and discuss how to derive, in 
the limit of infinitely large lattices, a finite set of coupled 
PDEs from the infinite system of ODEs of the previous 
section. Even if the PDEs are approximate, their solu- 
tions reproduce qualitatively the most probable (typical) 
behavior of the stochastic system at large times and large 
distances, the so-called hydrodynamic limit. In our case, 
this is not a very tight restriction since we have to dis- 
regard only a few time steps after the initial condition, 
and details on the scale of a single lattice step, such as 
the discontinuity close to the empty reservoir discussed 
earlier. On the other hand, the benefit is that only these 
PDEs will enable a discussion of the self-similar phenom- 
ena that takes place in the different models. 



A. Principle of the derivation 

We are interested in variations of the local cell density 
and local correlations on a length scale R much larger 
than the lattice step A. This length R can be the dis- 
tance between source and sink in a steady-state regime, 
or the size of the region of a Petri dish that cells exiting 
a spheroid have already invaded Fig. 2(b) If e M" 
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FIG. 6. Profile of |(a)| the density and |(b)| the connected cor- 
relation at different times for the adhesion model on a cylin- 
drical geometry with q=0.2 and L=75. The dots represent 
the results of the stochastic simulations and the error bars 
are smaller than the size of the symbols. The lines refer to 
the results of Eqs. ((9l 1121 I20|l . The insets show the scaling 
behaviors. 



In 



(b) 



the inset on the left-hand side refers to the 
analytical results, and that on the right shows the respective 
results obtained from cellular automata simulations. The pa- 
rameter = 6 is properly chosen to show the scaling of the 
stochastic simulations. 



is the position (in continuous space) of the center of lat- 
tice tile i, we assume that there exists some function p 
such that (rji) — p{ri/R) for all i and that p is differ- 
entiable as many times as necessary w.r.t. its argument 
Ti/R. This implies that the average number of cells in 
the lattice site i, (rji), varies from the average number in 
the nearest neighbor site j G V{i) from no more than a 
quantity of the order of 1/R. On the light of the previ- 
ous results, we know that it is an excellent approximation 
most of the time. The exceptions are in the set-up with 
migration out of a spheroid, at short times after the ini- 
tial condition when p has a steep slope, but this lasts less 
than sa 10 time steps, and in the steady state close to a 



FIG. 7. Profile of |7(a)| the density and |7(b)| the connected 
correlation at time t — 450 for the gap junctional model on a 
radial geometry with rare = 3 and Vsnk = 27. The lines show 
the results of 3 ■ 10^ averaged independent stochastic simu- 
lations at time t with same initial conditions and parameter 
p=l. The dots and the bars represent the average and the 
error for only 10^ stochastic simulations. 



sink of cells, for some particular values of the parameters 
of the models. 

Likewise, we assume that there exists a function C2 
such that the connected correlation function between 
nearest neighbor sites i and j e V{i) reads 



2R 



Since the correlation between the sites i and j does not 
depend on the order of i and j, it is natural to use the 
middle point of i and j as the first argument of C2 . The 
second argument of C2 can only take a finite number 
of values since the lattice is regular and it gives the di- 
rection of the two-point connected correlation function 
(sign is unimportant, i and j can be exchanged). To 
make expressions shorter, in the sequel, we will write 
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6*2 {r, Tj — Ti) as where / is the number of the direc- 
tion rj — Ti . Finally, we have seen in some of the previous 
simulations that finite regions of the system may tend to 
homogenize (the density p gets uniform). Therefore, we 
introduce the possibility that C' vanishes as R goes to 
infinity, and we assume that it involves a power law with 
the exponent ^ to be chosen later. The exponent ^ must 
be nonnegative since the correlation function is always 
between -1 and 1. 

It is generally not possible to define a regular func- 
tion of two arguments, say C2{ri/R,rj/R): indeed, 
C2{ri/ R,rji / R) may have a difference of order 1 (not 
1/i?), if j' is another nearest neighbor of i than j, but 



without rji — ri being colinear to Vj 



This is for 



instance the case on the cylinder of Fig. [5] where and 
(e.g. correlations along the cylinder axis and perpen- 
dicular to the cylinder axis) can be quite different one 
from another because the boundary conditions break the 
rotation invariance of the lattice. The same holds if f is 
not a nearest neighbor of i. 

Then we insert the expressions for p and C2 defined 
above into the set of discrete equations p4)) for the time 
evolution of p{i) and C2(«,j), and we perform a Taylor 
expansion in powers of R~^ on the right-hand side of the 
equations around the points ^ and ^^2R^ respectively, 
neglecting terms of orders three and above in R^^ . There 
are 1 + d equations for each lattice site, where d is the 
number of different directions, hence a total of (1 -I- d)N 
equations for the whole lattice of N sites. But, since the 
equations are independent of the site up to an index shift, 
except for the sites on the boundaries of the system, we 
end up with only \ -\- d PDF, plus boundary conditions 
— the (approximate) "continuous space model" , or "hy- 
drodynamic limit". 

Of course, one could systematically generalize this pro- 
cedure to take into account correlations at distances 
larger than one lattice step or between three or more 
points — this may further improve the quality of the ap- 
proximation of the continuous space model. To do this, 
one can redefine the multi-points density function: 



Pp(f7(n), • • ■ , r]{rp)) = pp(7?(Ari), . . . , r?(Arp); rc) 



where rc is the centroid — , and Ar^ = rk — rc 
is the position of k-th point of pp with respect to the 
centroid. Under the assumptions above, pp with fixed 
Arfc for all fc is a regular function of rc- Then we could 
perform a series expansion, at fixed Ar^, in powers of 
1/i?, of the multipoint probability density which can be 
in turn inserted into the system of discrete equations to 
yield a system of PDF. But here, we restrict ourselves to 
equations with p and only. 



B. Explicit expressions 

Before explicitly applying the Taylor expansion to the 
discrete equations for the evolution of the density 
and the connected correlations, it is important to remark 
that microscopic symmetries related to the regularity of 
the lattice reflect on the hydrodynamic limit. 

Therefore, the equation of the two-point correlation 
between neighbor sites in one direction can be changed 
in the equation for the correlation in another direction 
by applying rotations of 7r/3 rad as consequence of the 
invariance of the hexagonal lattice under such transfor- 
mations. This symmetry holds everywhere in the region 
of interest between the reservoirs, except in the sources 
and in the sinks. 

To easily express the results and the symmetries of the 
system, we introduce a set D of unitary vectors identify- 
ing the three main lattice directions eo, ei, 62, S such 
that: 



eo 




62 




(23) 



and the directional derivative = e^.V, where the index 
i means that the derivative is taken along the direction 

In the case of the linear model, after performing the 
expansion in series in respect to the lattice distance A to 
each term of the right-hand side of Fq. (HH) , the results 



dtp{r,t) = 



B. 



k,p 



i?2 



3 

E 



B 



k,Ci 



B 



k,p,Ci 



, (24) 



dtC^{r,t) 



fee{a,/3,7} . 



Ri ' i?2 



3 

-E 

1=1 



B^ 



B^ 

k,p,C^ 



(25) 

where A is the zero order and the terms named B are the 
second order of the Taylor expansion. In both equations, 
for the density and the correlations, there are no terms 
with odd powers of 1/i? due to reflection symmetries of 
the lattice with respect to the lattice main directions. 
The equation for the density does not have any zero or- 
der term because the number of cells is locally conserved. 
The terms i3fc,p, and B^- c^ both contain terms propor- 
tional to the parameter k e {a, /3, 7} and each term only 
depends on the density and the correlation in direction 
i, respectively. Bf^ p c^ is the interaction part that takes 
into account the coupling between the density and the 
correlation in direction i. The same definitions hold for 
the terms in Fq. (|25p where the upper index i (resp. j) 
refers to one of the main directions of the two-point con- 
nected correlation. From here on, all upper indices must 
be considered modulus three, even though it is not explic- 
itly written. Using the notation previously introduced, 
the A and B terms in Fqs. ([Ml [SS]) are: 
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B. 



= B. 



B. 



Ba,p,Ci + ^ — Ba,p,Ci+^ = 



Bb 



B^,co^j^{^A-4Vt)C' 

i?^,ci = ^(3A-4V?)Ci 

Bp,c^ = ^{3A-4S7l)C^ 
B0,p,c° = Bp^p^c^ = Bp^p^c^ = 



^ j2(5 - 7p)[(Vip)2 + V2PV1P + (V2P)'] + ^(10 - 7p)pAp 



B 



7, CO 



7,Ci 



B 







^7,P,co = [(1 - P)(4Vo + 3A)C° + C°(3A - 14Vg)p - 18VopVoC°] 



B^,p,c^ 



12 
1 

12 



[(1 - p)(4Vf + 3A)Ci + C\3A - 14Vf)p - ISVipViC^] 

[(1 - /9)(4Vi + 3A)C2 + C2(3A - 14V^)p - I8V2PV2C2] 



Al = -(C^+i +C-'+2 -5C^) 

O 

-i[(V,+l+V,+2)p]' 



a,p 



= B-' = f) 



a,C5 + i ~ 24 ^'+2 



= -[C^+i + - 2(1 + 3p)&] 
= ^ [(1 - 4P)(V,P)' - 2p(l - p)V,+iV,+2p] 



B^ 



Bi 



B^ 

I3,P,C3 



6 




24 



'i+i 



(0,+2 



^C-'(28V,+iV,+2 



57A)p 



Bip,c^+i = ^[C'+''^j+i - (V,+2C^+^)](V,- + V,+i)p 
^1p,c.+= = ^[^^'■+'V,-+2 - (V,-+iC^+2)](V,- + V,+2)p 
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Ai> = 3 {P [-3p(l - pfR-^ + {P + 2)(C^+i + - (7p + 5)C^] } - 

Bip = ^ { (2V - 26p + 3) {V,pf + p{p - 1) [(3 - 14p)(V,Vi + ^%^) + 6(1 - 4p)(V,+iV,+2)] p] 

^;,P,c. = ^ {4(1 " P)PV,2C^ + 4(4p - 3)V,C^V,p 

-0 [l46(Vjp)2 - 224Vj+ipVj+2P + (21Vj2 - 32Vj+i Vj+2)p + p(158Vj2 - 128Vj+i Vj+2)p] } 
^7,P,c.+i = ^ {P(P- 1)(V,\2 - 8V,+iV,+2 - ^V]+i)C^+' + 3V2+2C-'"+i + 2C-'"+iV,p(V, + 16V,+2)p 
+2(1 + p)C^+i(26Vj+iVj+2 + SV^^i + 9V2^2)P - 64C-''+iVj+iVj+2p 
+8Vj+iC^+i[(2p - 1)V, + (p - l)V,+2]p - 4Vj+2CJ'+M(p + l)Vj + 4(2p - l)Vj+2]p} 

^1p.c.+^ = ^ {P(P- l)(V,\i - 8V,+2V,+i - 8V2+2)C^^+' + 3V2+iC^"+2 + 2C^"+2v,p(V, + 16V,+i)p 
+2(1 + p)C^+2(26Vj+2V,+i + 5V2+2 + 9V2+i)P - 64CJ'+2Vj+2Vj+ip 
+8Vj+2C^+2[(2p - 1)V, + (p - l)V,+i]p - 4Vj+iCJ'+2[(p + l)Vj + 4(2p - l)V,+i]p} 



Defining S and s the spatial subsets of M? where there 
are sources and sinks respectively, the border and the 
initial conditions used to determine the solution of the 
system of partial differential equations are: 

p(r,t) = l Vre5, 

p{r,t) = Vr e s, 

C*(r,t) = VresUS', 

p(r, i = 0) = C*(r, < = 0) = Vr. 

The computation of the hydrodynamic limit for the 
adhesion model and the gap junctional model do not 
present ulterior problems, but a subtle more analytical 
complexity, and it can be performed in the same way as 
for the linear model. The hydrodynamic limit introduces 
some approximations in comparison to the discrete equa- 
tions because one disregards terms of order bigger than 
1 /R^ in the series expansion. The numerical solutions of 
the PDEs are obtained by first re-discretizing the space 
derivatives using a second order finite difference method 
and then following the same numerical procedures as for 
the ODEs. They are accurately in agreement with the 
solutions of the ODEs for all the models proposed except 
when close to the sink. The analytical solutions are 
therefore, at the last step and two lattice steps away 
from the empty reservoirs for the density and the cor- 
relations respectively, the numerical solutions smoothly 
go to zero in all cases making it impossible to retrieve 
the discontinuity observed in the discrete representation 
for some values of the parameters of the models at the 
steady state. 



C. Self-similar behaviors 

When non-interacting agents randomly jump to the 
nearest neighbor site with a constant transition rate, in 
the continuous limit, the concentration p of agents at 
the position r obeys: dtp{r,t) — DAp{r,t). A gener- 
alization of the previous equation is the diffusion equa- 
tion in porous media [tI]: dtp{r,t) — W[D{p)'V p{r,t)] 
which is useful to describe the concentration of agents 
when the agents interact [tII- Both equations with the 
initial condition p{r,t = 0) — 5{r) show a self-similar 
behavior. Indeed, the solution of the diffusion equa- 
tion at a given time ti is self-similar to the solution at 

time 12'- p{r,ti) — p{r ^J^,t2)- In this specific case, 

the self-similarity of the density function is exact, but in 
other cases, when correlations are taken into considera- 
tion, such behavior holds true only approximately and at 
large times far from the initial transient regime [75| . For 
the geometrical dispositions of the reservoirs proposed in 
this work, the presence of a sink destroys the self-similar 
behaviors of the density of cells; therefore, we consider 
the situation where the sink is sufficiently far (T < L, 
where T is the maximum time we run the stochastic sim- 
ulation) so that the cells can migrate for long enough 
without perceiving the sink. 

To investigate the self-similarity of both the density 
and the two-point connected correlation, we exploit the 
PDEs derived from the hydrodynamic limit for the lin- 
ear model to retrieve the scaling of the solutions. When 
there are no correlations, the Eq. ((M)) becomes a particu- 
lar type of diffusion equation in porous media; therefore. 
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after eliminating the factor thanks to the change of 
variable t = i, we can write the solution of Eq. ([M)) 
as: 



r 



(26) 



where the function g{r,t) on the right-hand side tells 
us how the exact solution of the density profile differs 
from the self-similar result /("^) obtained in case of no 
correlations. The parameter e is the magnitude of such 
discrepancy and for negligible correlations in comparison 
with the density, it holds true that e ^ 1. Inserting the 
solution /(-^) into Eq. (pS)) . one can find the approxi- 
mate solutions for the correlations and use them again 
in the equations for the density. Repeating the process 
iteratively, more accurate solutions for the density and 
the correlations can be found, and at each repetition, it 
is possible to check the validity of the hypothesis e ^ 1. 
In the insets of Figs. 3(a)[ 4(a)[ 6(a) the stochastic sim- 
ulations of each model at different times are scaled using 
the scaling behaviors of the function f{-^)- We have 
observed that the same self-similar behavior holds for 
stochastic simulations of the linear model with various 
values of parameters which reproduce the adhesion or gap 
junctional model behaviors. The perfect overlap between 
the curves shows that analyzing the scaling of p{r, t) un- 
der the conditions of negligible correlations is very good. 

In the sequel, we change the notation t to t for sim- 
plicity, and we consider that p{r,t) — /(-^)- Eq. (l25l) 
now reads: 



fce{a,/3,7} . 



Ai 



3 



i=l 



k,p,C' 



(27) 

Linear model — only exclusion interactions. To 

begin, let us consider the simplest case a ^ and (3 = 
7 = 0. In Eq. (P7|) . choosing ^ 7^ 2 leads either to a 
trivial solution C = or to an inconsistent solution V/o = 
0. Therefore, the only possible choice is ^ = 2 which 

results in C^{r, t) ^ Ap{r, t) 



while 9tCJ(r,t) 



S/'jC^r,t) are asymptotically negligible. 

Linear model — mimicking gap junctions p < 1. 

The same discussion is valid for 7 = and both a and 
/3 different from zero. This brings us to the same self- 



similar behavior (r, t) ^ Ap{r,t) 



In the left 



insets of Fig. |4(b)[ the numerical solutions of Cy for the 

gap junctional model are rescaled using C-' (r, t) — ^ 
as self-similar behaviors. The good agreement given by 
the overlapping of the curves also shows that the gap 
junctional model for p < 1 share the same self-similar of 
the linear model. 

The right inset of Fig. |4(b)| shows that the results of 

the stochastic simulations scale as C^{r,t) ~ ln(t) — 
The logarithmic correction in the scaling of the correla- 
tions is due to the long range correlation produced by the 



exclusion process, and they cannot be retrieved using the 
hydrodynamic limit equations. 

Linear model — mimicking gap junctions p = 1. 

The case a = 7 = and /3 7^ is particularly interesting 
because the solution of Eq. ([M)) with no correlations is 
no longer a smooth function of the position, but it is a 
piecewise function with two regions: the part closer to the 
source, where the density decreases with a finite negative 
slope as its distance from the full reservoir increases, and 
the part that goes from the positions where the density 
becomes zero onwards, in which the density is uniformly 
null [7^. 

The discontinuity in the first derivative of the density 
profile for a = is true only when there are no cor- 
relations, while, when correlations are introduced, such 
discontinuity disappears. Hence, always under the con- 
dition p{r,t) = /("^): in the region where p{r,t) ^ 0, 



the correlation rescale as C^{r,t) = 

On the other hand, in the region where the density is 
zero, all terms containing p or its derivatives are zero, 
and Eq. (P7|) becomes: 



/3 



4 



3 

-E 



(28) 



The leading order term in 1/i? must be zero. Hence, 
A^p = Ap = — and Eq. (^5)) becomes a system 
of three standard diffusion equations without sinks or 
sources. From the expression of A^, we conclude that 
C° = = C^. Therefore, Eq. ^ simplifies to a stan- 
dard diffusion equation for any of the quantities . If 
we look for a nontrivial result for , the three diffusion 
equations become identical to: dtC^ — ^AC'^ . 

From the conservation of the correlation and the 
solution of the diffusion equation follow that the correla- 



tions in all directions scale as: C^(r,t) 



If the 



complete solution Eq. ([^ is used instead of the approx- 
imate solution of the density with no correlations, then 
the region with p = changes in p <C 1. This correction, 
even though it is very small (e <^ 1), solves the incon- 
sistency of having a correlation different from zero in a 
region with no cells, and leaves unchanged the self-similar 
behavior. 

In the top and bottom left insets of Fig. |3(b)[ we see 
that the numerical solutions of the ODEs at different 
times for the gap junctional model with p = 1 reproduce 
the self-similarity of the case a = 7 = 0, and (3^0, 

where the negative tail of the correlations scales as — — 

and the positive peak of the correlations as jl* 
top and bottom right insets of Fig. |3(b)| show the scaling 
for the stochastic results of the gap junctional model with 
p = 1, and also in this case the logarithmic correction is 
due to the long range correlations. 

Linear model — mimicking adhesion. In the case 
with the parameter 7 7^ 0, the scaling behavior is given 



t 

The 
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by which contains terms proportional to positive pow- 
ers of p{r,t) and R^. One has to choose ^ = to get 
consistent results. Consequently, the self-similar form of 
the correlation is C^{r,t) ^ p{r,t) ^ f(^)' meaning 
that the height of the correlation peak remains constant 
during the migration process. When both /3 and 7 are 
negative, the linear model mimics the behaviors of the 
adhesion model. Indeed, scaling the numerical solutions 
at different times of the ODEs for the adhesion model 
with the same scaling behavior of the corresponding lin- 
ear model, we have a perfect overlapping of the curves, 
left inset of Fig. [6(b)] In the right inset, the results of the 
stochastic simulations show an increasing of the height of 
the peak due to long range correlations. In the adhesion 
model such effects are much smaller and asymptotically 
negligible than in the other cases where correlations pro- 
duce logarithmic correction. In all the insets in Figs. [31 SI 
[6l the curves that deviate the most from the self-similar 
behaviors are those where the sink has been reached and 
the constraint t < L does not hold any more. 



V. DISCUSSION 

A. Measuring the correlations in experiments 

Let us try to estimate for what parameters (popula- 
tion size, number of repetitions, duration of experiment) 
one could get an exploitable measure of the connected 
correlation function in migration assays of cells on Petri 
dishes. To be closer to this experimental situation, we 
performed stochastic simulations and got numeric solu- 



tions of Eqs. PHITTI) on the radial geometry of Fig. 2(b) 



In this geometry, the density profile is computed by av- 
eraging the densities over the sites belonging to the same 
annulus centered on the source center. Ogre, and one lat- 
tice step thick, so that the density depends only on the 
distance from the source center. The profile of the near- 
est neighbor two-point connected correlation is obtained 
by averaging the correlation over the six nearest neigh- 
bors of one of its sites, then the result is averaged over 
all the sites belonging to the respective annulus as for 
the density profile. At large times, the radius of the 
outer region invaded by cells is large, so that one can 
expect that values of the observables in that region will 
be similar to what is found with the cylindrical geome- 
try of Fig. 2(a) It may be different both at short times 
and close to the center of the disk. However, it turns 
out that the "correlation wave" observed and discussed 
in Sec. pil[) and Sec. (jIVp is still present in this setup 
and even amplified. In addition, in the case of the gap 
junctional model where correlations are the most impor- 
tant far from the center of the disk, the relative error on 
the measure of the correlations tends to decrease with 
time, since the measure is done by averaging over a re- 
gion which gets larger and larger. To show these results, 
in Fig. [71 we have plotted the average of 1000 repeated 
simulations at a given time and the respective error bars 



in comparison with the results obtained from the aver- 
age of a much higher number of simulations for the gap 
junctional model. The figure, on one side, shows that it 
is not necessary to have a prohibitively high statistic to 
observe the two-point connected correlations, while, on 
the other side, it shows that simulations on the radial ge- 
ometry are very useful to compare the proposed models 
with experimental data on Petri dishes. 



B. Perspectives 

Let us comment on the connection of the present ap- 
proach to other techniques. The structure of Eq. (j26|) can 
be seen as an iterative computation of the local density p, 
where the first step would yield the mean-field approxi- 
mation to the evolution equation of p (all correlations C2 
being replaced with 0) and the second step involves an 
improvement to this equation that is expressed by means 
of C2 . If one formally integrates the evolution equation 
for C2 , replacing in it the value of p that solves the mean- 
field equation for p, C-2{t) gets expressed as an integral 
of p{t) over the time range < t' < <. Then, replacing 
C2(t) in the equation for p yields an integro-differential 
equation for p[t) — a non-markovian model with mem- 
ory kernels [77] . This is much like what is obtained in the 
Mori-Zwanzig formalism [TSj , the non-markovian charac- 
ter being due to a tentative of description of a complex 
situation (many spatially inhomogeneous configurations) 
with a simple quantity (the local density p) by a kind of 
"projection" . However, there are many ways to do this 
systematic reduction of numbers of freedom [U [77l - [79j , 
and, when truncated to the first or second step of itera- 
tion, not all schemes yield analytical approximation with 
the same degree of agreement to simulations. 

It would also be interesting to have a deeper under- 
standing of which situations (and possibly a criterion to 
discern them) allow one to get analytical results of good 
quality while taking into account only short-range cor- 
relations. Our situation might be related to the case of 
the Smoluchowski theory of aggregation discussed in [2l[ . 
The success of this method in spite of the simple ap- 
proximation it relies on may be explained, at least for 
some of the models it was applied to, by the observation 
that there is no propagator renormalization in the corre- 
sponding field theories (in an RG approach), and hence 
no anomalous dimension for the diffusion constant or the 
fields. 

There are natural and biologically relevant extensions 
of the models considered in this work that can be stud- 
ied with the same techniques, and we plan to deal with 
some of them in our future work. First, one can incorpo- 
rate the effect of cell proliferation. According to previous 
works [sll - fsH , cell proliferation leads to short range cor- 
relations (because daughter cells are close to the position 
of the mother cell) which have the same order of mag- 
nitude as the ones produced by our contact interactions. 
In a realistic model, one should probably take care of 
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the migration-proliferation dichotomy in the behavior of 
cells [io, HH, which may enhance correlations. Then, the 
technique should probably be useful for other models of 
exclusion processes with interactions like the ones of [82j 
and [HI. It can also straightforwardly be extended to 
3D migration, for which experimental data is more diffi- 
cult to obtain but still accessible thanks e.g. to confocal 
microscopy. 

An extension to disordered systems would also be use- 
ful, for instance in the case of deformed lattices [H, [11] , 
since they yield more realistic models of biological pro- 
cesses and make it possible to avoid some regular lattice 
artifacts. 



C. Conclusion 

We have shown how to extend simple mean-field ana- 
lytical approximations of spatially inhomogeneous exclu- 
sion processes with local contact-like interactions as can 
be found in biological situations of interacting migrat- 
ing cells. Supplementing the local density of agents (or 
cells) with short-range correlations in the analytical, de- 



terministic, macroscopic approximations of the stochas- 
tic processes yields not only results in better agreement 
with stochastic simulations. It can also be used to in- 
fer more precise estimates of the interaction parameters 
from experimental data about density and correlations 
of migrating cells than from density values alone, when 
the type of cell-cell contact interactions is known. If the 
latter is unknown, it can be discovered by comparison 
with the classification scheme that we provide, which was 
made possible thanks to the PDEs approximation of the 
exclusion processes: the self-similar behavior of the two- 
point connected correlation function is highly dependent 
on the type of contact interactions, adhesive or contact- 
maintaining for instance, at the microscopic level. 
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